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ARTICLE INFO ABSTRACT 


Keywords: Hypothesis: The large-scale implementation of hydrogen economy requires immense storage spaces to facilitate the 
DLVO Theory periodic storage/production cycles. Extensive modelling of hydrogen transport in porous media is required to 
Wetting comprehend the hydrogen-induced complexities prior to storage to avoid energy loss. Wettability of hydrogen- 
i a brine-rock systems influence flow properties (e.g. capillary pressure and relative permeability curves) and the 


residual saturations, which are all essential for subsurface hydrogen systems. 

Model: This study aims to understand which parameters critically control the contact angle for hydrogen-brine- 
rock systems using the surface force analysis following the DLVO theory and sensitivity analysis. Furthermore, 
the effect of roughness is studied using the Cassie-Baxter model. 

Findings: Our results reveal no considerable difference between Hz and other gases such as Np. Besides, the in- 
clusion of roughness highly affects the observed apparent contact angles, and even lead to water-repelling fea- 
tures. It was observed that contact angle does not vary significantly with variations of surface charge and density 
at high salinity, which is representative for reservoir conditions. Based on the analysis, it is speculated that the 
influence of roughness on contact angle becomes significant at low water saturation (i.e. high capillary pressure). 


1. Introduction 
1.1. Underground hydrogen storage 


Fossil fuels have been the primary energy source since the beginning 
of the industrial revolution due to their large (volumetric) energy den- 
sity, accessibility, ease of production and transport. However, they are 
also the main cause of climate change, because of their large carbon 
footprint. The transition towards low-carbon renewable energy sources is 
inevitable to reverse global temperature rise and to avoid further harm to 
the Earth. A wide array of alternative energy resources are available 
causing very low carbon emissions. However, their major drawback is 
their intermittent nature that gives rise to periods of energy shortage and 
surplus throughout the year [1]. The excess electricity produced by these 
renewable sources during a productive season could be stored for reuse at 
times of need. The excess energy can be stored at the surface, for 
example, in high-pressure gas cylinders, cryogenic tanks, and batteries. 
The South Yorkshire battery storage and the upcoming Greater Man- 
chester liquid air storage sites are good examples of such storage systems, 
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holding about 300 MWh of energy [2]. However, to meet the (ever-- 
increasing) global energy demand, storage in the order of TWh is 
required. 

From “Power to Gas to Power” framework, the excess produced 
electricity can be stored in the form of hydrogen gas [3]. Hydrogen can be 
used as an energy carrier as it has a high energy density (120 MJ/kg) and 
is easily convertible to other energy forms. However, this scheme re- 
quires large storage spaces to hold hydrogen for different time intervals 
(days, months, or even years). Subsurface formations may provide a so- 
lution as the porous underground rocks (aquifers, and depleted gas res- 
ervoirs) can facilitate the safe and economical storage of hydrogen gas. 
The research on the underground gas storage (e.g., natural gas, CO2, 
compressed air) is rich [4-7], nevertheless, hydrogen exhibits unique 
physical and chemical characteristics compared to other gases that can 
cause various operational problems and lead to energy loss reducing the 
overall efficiency of the process. The inherent properties include low 
density, low viscosity/high mobility, and high reactivity. Consequently, a 
comprehensive view of the hydrogen-induced processes in the subsurface 
formations is required prior to large-scale field applications. 
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1.2. Wettability: surface versus porous media 


Wettability is one of the material properties which controls flow and 
transport in porous materials. The conventional way of quantifying 
average macro-scale wetting behaviour of a three-phase system has been 
through contact angle measurements/calculations which usually take 
place on a flat surface using methods such as sessile drop and captive 
bubble [8]. The upscaled effect of wettability can be estimated using 
capillary pressure-saturation curves by methods such as USBM and Amott 
tests [9,10]. More recently, researchers have focused on predicting the 
wetting properties of porous systems through X-ray and micro-CT im- 
aging at different scales. In some of these studies, the fluid surface 
coverage parameter determines the wetting state of the rock rather than 
the three-phase contact angle [11-13]. The key feature of these studies is 
considering the two-phase flow conditions to derive the in-situ wetting 
state. Another newly-developed methodology uses the Gauss-Bonnet 
theorem and deficit curvature definition to derive wetting [14-17]. 

In the majority of studies related to wettability, highly-polished fine 
surfaces are used. However, surface roughness has a considerable impact 
on the measured/calculated contact angles and on contact angle hys- 
teresis [18]. For instance, on a hydrophilic surface, the increased surface 
roughness leads to a more water-wet system, while for hydrophobic 
surfaces, roughness increases the chance of hydrophobicity [19,20]. 
Although the majority of the previous studies reported a similar trend 
[21-28], there exists a few studies challenging it [29]. Kaveh et al. [30] 
studied the wettability of a Bentheimer sandstone at the presence of CO2 
and flue gas (CO2+Ng). Their results demonstrate a wider contact angle 
range on the less-polished surfaces and higher chance of hysteresis. The 
experimental results of Hosseini et al. [22] also revealed inverse rela- 
tionship between the contact angle and calcite surface roughness. Sari 
et al. [31] concluded that the increase of pore roughness reduced the 
effectiveness of low-salinity water (LSW) injection in modifying the 
wettability of carbonates [31]. Alnoush et al. [21] conducted AFM 
analysis to investigate the contact angle variations on the 
naturally-cleaved and mechanical-ground calcite surfaces and concluded 
that reducing the diamond grit size (increased roughness) leads to 
reduced water contact angles (more water-wet) [21]. Tudek et al. [29] 
performed a wettability study combining sessile drop and CT-Imaging 
techniques along with Lattice-Boltzmann modeling for liquid CO2-bri- 
ne-sandstone systems. As opposed to the previous studies, their 
sessile-drop contact angles revealed a direct relation to the surface 
roughness with an average of 30°. 

The contact angle of the hydrogen-brine-rock systems was measured 
using captive bubble [32,33], tilted plate [22,34], and other methods 
[35-37] as summarized in Table 1. However, results do not provide a 
conclusive trend or value and it raises the question whether discrepancy 
in the results is due to multiple factors being different across the exper- 
iments Muhammed et al. [38]. 

From the rich literature on wettability and recent advances in this 
field, four different factors can be counter responsible for the contact 
angle in porous media: surface forces due to chemistry of the system, 
surface roughness, three-dimensional pore geometries, and dynamic 
conditions. In the present study, we will focus on the first two groups 
(surface forces and surface roughness) to identify theoretically what are 
the critical physical-chemical parameters. 


Table 1 
The reported literature experimental data on hydrogen-brine-rock 
contact angles. 


Experimental Study Contact Angle Range 


Yekta et al. [35] 21-35° 
Hashemi et al. [32, 33] 21-43° 
Al-Yaseri et al. [36] 10-40° 


Ali et al. [34] 5-65° 
Hosseini et al. [22] 30-80° 
Iglauer et al. [37] 5-50° 
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1.3. The present study 


This study aims at evaluating the hydrogen-brine-rock wettability 
relevant for underground hydrogen storage. We particularly focus on 
investigating whether hydrogen poses additional complexities compared 
to other gases, such as N2, CO2, and CH4. Hence, we utilize the surface 
force analysis to investigate the impact of different parameters on the 
rock wettability in the presence of hydrogen and water. Although a 
similar analysis has been performed for petroleum-related problems, the 
inconsistency in the literature data justifies the application of this 
method to describe the wettability of hydrogen-brine-rock system, 
especially when small-scale surface roughness is included. This study 
provides the following key contributions by application of DLVO theory 
and surface roughness analysis: 


@ Probing into how main physical and electrochemical properties of the 
rock/water affect the contact angle on smooth surfaces. 

@ Evaluating the influence of subsurface thermophysical conditions on 
the contact angle. 

@ Exploring how the surface roughness affects the apparent contact 
angles during the drainage process. 

@ Investigating whether hydrogen imposes unique complexities to the 
storage process compared to other gases, like No. 


2. Surface force analysis and contact angle calculation 


The literature shows contradictory results on the impact of parame- 
ters such as pressure and temperature on the contact angle. Furthermore, 
experimental studies on wettability are challenging, especially for highly- 
flammable gases like hydrogen or supercritical fluids like CO2. Given the 
uncertainties and contradictory results, there is a need for robust theo- 
retical models that can help with understanding the experiments. As 
such, wettability and the resulting contact angle can be quantified by 
means of analytical surface force calculations [39]. 

There are three main force types affecting the particles in a colloidal 
dispersion or solution. The van der Waals (vdW) force is the first force 
type acting on the particles and by definition is always attractive between 
similar objects. The electrostatic forces, on the other hand, are respon- 
sible for retaining these particles in the solution by creating a repulsive 
electric field. The interaction of these two forces forms the basis of the so- 
called DLVO (named after Boris Derjaguin, Lev Landau, Evert Verwey, 
and Theodoor Overbeek) theory [40]. It has been shown that the con- 
tinuum theories of the DLVO forces are only applicable to particle sep- 
arations above the molecular diameters. When the separation falls below 
a few nanometers, these theories fail owing to the presence of the 
non-DLVO structural or hydration forces that will come into play at short 
distances. Hydration forces will be an addition to the repulsive in- 
teractions for hydrophilic particles and are caused by layering of water 
molecules between the two surfaces [41]. One of the first studies on 
estimation of wettability using DLVO theory is the work of Hall et al. 
[42], in which they performed an analysis of wetting film stability on 
Athabasca Oil Sands and Buckley et al. [43] who studied the adhesion 
properties of several crude oils on solid glass surfaces. The DLVO theory 
was further used in other subsurface applications such as thin film sta- 
bility [44], calcite wettability [45], wettability during low-salinity 
waterflooding [46-53]. 

Assuming a triple layer system (hydrogen-brine-rock) on a flat surface 
at the thermodynamic equilibrium conditions (Fig. 2), the capillary 
pressure should be at balance with the disjoining pressure in the thin 
water film, i.e., 


P. = I(ho), @) 


in which, P, is the capillary pressure, II denotes the total inter- 
molecular force per unit area between two surfaces (also referred to as 
disjoining pressure) and họ refers to the equilibrium film thickness. The 


F. Nazari et al. 


e van der Waals 


Disjoining Pressure 
© 


+ Electrostatic 


= = Structural 
Film Thickness 


Fig. 1. Schematic presentation of three components of total disjoining pressure. 
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Fig. 2. A schematic of the film/meniscus transition region adapted from Hir- 
asaki [39]. At film thicknesses larger than the equilibrium one (ho), the contact 
angle changes with film thickness until it becomes constant at equilibrium 
contact angle. 


disjoining pressure is also the sum of the van der Waals, electrostatic 
double layer, and hydration contributions as demonstrated in Equation 


(2). 


Tota’ = Maw + Hent + Hyr (2) 


Note that negative values indicate attractive interactions, and positive 
ones indicate repulsion. A schematic of the individual components of the 
total disjoining pressure has been shown in Fig. 1. The first section in the 
supplementary material file presents the detailed methodology for 
calculating different components of the total disjoining pressure based on 
the DLVO theory. 

A schematic of the film and meniscus region is shown in Fig. 2. For 
any film thickness larger than then equilibrium film thickness, contact 
angle, @m would change with film thickness following cos ôm = 1— 


| i (P: —I1)dh| 1 o, where o is interfacial tension (IFT). The equilibrium 
contact angle (6¢,) can also be calculated as follows [39,54,55]: 


cos 0y = 1 — | [nn + hat) /s (3) 
ho 


The vdW, electrostatic, and hydration contributions were calculated 
using the formulations described in the supplementary material file. The 
Lifshitz theory was employed for calculating the Hamaker constant. 
Using the calculated disjoining pressure, equilibrium contact angle (¢q) 
was determined from Equation (3). The equilibrium film thickness (ho) 
was found by intersecting the total disjoining and capillary pressure 
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(assigned based on the reservoir conditions). This led to a nonlinear 
equation which was solved using the Newton-Raphson method. Note that 
the contact angle calculated from Equation (3) is driven for smooth 
surfaces. The contact angles on a rough surface (8a) can be several times 
higher. To account for the presence of gaseous fluid and surface rough- 
ness, the Cassie-Baxter model can be implemented [19,56—59]: 


cos 0, = r;(ficos 0; + focos 62) (4) 


In this equation, 64 is the apparent/rough surface contact angle, rp is 
the roughness ratio, f; and f are fractions of solid surface wet by the 
water and hydrogen, respectively [19,20]. In order to determine f and fo, 
a simple V-groove geometry as demonstrated in Fig. 3 has been consid- 
ered. The protrusion depth (d) and spacing (w) between its peaks can be 
determined by Atomic Force Microscopy (AFM) measurements [22,31, 
60]. Using these parameters, one can calculate the fraction of the surface 
that has been wet by the water and hydrogen as a function of the inter- 
face curvature (r) and capillary pressure. Note that the wet water area 
shrinks as the capillary pressure rises (water saturation decreases) within 
the system. From this point forward, the smooth surface contact angle 
derived from the DLVO calculations will be referred to as “equilibrium 
contact angle” and the one on rough surfaces will be referred to as 
“apparent contact angle”. 


3. Results and discussion 


First, in order to ascertain the reliability of the calculated film 
thicknesses and other dependant variables, we have benchmarked our 
model against the paper of Radoev et al. [61]. The present model leads to 
an equilibrium film thickness of about 22 nm at 0.01 M salinity, similar to 
the results presented at Fig. 7 of Radoev et al. [61]. After checking the 
validity of the model against the reported data, the hydrogen-brine-rock 
system was analysed. 

The initial parameters for the model are shown in Table 2. The sub- 
scripts 1 and 2 denote the rock-water (RW) and water-hydrogen (WH) 
interfaces, respectively. The potential for the hydrogen-water interface is 
taken after the reported values of Takahashi [62] for air-water systems at 
around neutral PH values. The surface charge density for this interface is 
also assumed due to the lack of data. The refractive index and dielectric 
constant of the rock surface are also selected based on the Butt et al. [8] 


_ 2h + fw) 


fo w 


Fig. 3. A schematic of the assumed surface protrusion and its corresponding 
parameters. The ry, fi, f2, and r will be calculated as a function of w, d, and Beg 
based on the presented equations. 
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Table 2 

Default parameter values used in this study. 
Parameter Value Parameter Value 
€n (-) 4.50 [8] 61, C/m? —0.004 [63] 
n, (—) 1.50 [41] 02, C/M? —0.004* 
nw, (—) 1.33 [41] Ve, Hz 3.00 x 1015 [8] 
Mp, (—) 1.00 [64] Temperature,K 340.00 
Wi, mV —10.00 [65] IFT, mN/m 60.00 [66] 
Wo, mV —20.00 [62, 67] Pa kPa 110.00 [35] 
Salinity, M 0.10 Salt MW 58.44 
Roughness (—) 1.33 [29] 
* Assumed 


for quartz surfaces. 

Fig. 4 shows an example of the calculated disjoining pressure iso- 
therms and the local meniscus angle for a hydrogen-brine-rock system 
assuming constant potential boundaries. The film thickness and equi- 
librium contact angle for the case illustrated in Fig. 4 are about 2.14 nm 
and 6.22°, respectively. It is also apparent that at very small separations, 
the electrostatic contribution is negative and becomes positive at larger 
separations. This is consistent with the observations of Hirasaki [39]. The 
reason is that when the surface potentials have the same sign, the surface 
with a higher potential value provides a higher contribution to the 
electrical field between two bodies, and the second surface goes through 
a charge reversal [39]. 


3.1. Sensitivity analysis 


In this section, a sensitivity analysis is performed to evaluate the ef- 
fect of different model parameters on the wetting features of the 
hydrogen-brine-rock system. The sensitivity variables are the surface 
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Fig. 4. (a) The disjoining pressure isotherms and (b) the meniscus angle of a 
sample hydrogen-brine-rock system. Results are for a constant potential case 
with the model parameters shown in Table 2. As it can be seen, the vdW 
isotherm dictates the total pressure and the type of the interaction. 
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potential, surface charge density, capillary pressure, interfacial tension 
(IFT), temperature, salinity, and dielectric properties. Where applicable, 
the plots will be generated at two largely-varying salinity values (0.01 


and 1.00 M NaCl). The range of parameters used in the sensitivity study 
are presented in Table 3. 


3.2. Surface potential 


Fig. 5 illustrates the outcome of the conducted sensitivity by varying 
the surface potential for the two salinities. The potential range is 
considered such that it can encompass both WH and RW values (here 
from —60 to 10 mV). The second interface potential is set to an average 
value of —15 mV [43,62]. As can be seen from the graphs, increasing the 
salinity reduces the magnitude of the electrostatic pressure and therefore, 
the meniscus angle variations are much lower for the high-salinity case. 
This is attributed to the elevated Debye length caused by increased ionic 
strength [8,41]. 

At low ionic strength, when the surface potentials have the same sign 
(—40 to 0 mV), the electrostatic isotherm is negative at small separations 
and it becomes positive at larger separations. When the surface potential 
changes sign (0-10 mV), the isotherm becomes negative/attractive at all 
separations [39]. However, the case with high ionic strength or salinity 
behaves differently, namely, 1) the negative isotherm starts from much 
lower separations (Fig. 5b inset), 2) for the average potentials (mainly 
around —15 mV), the isotherm is positive across the board, 3) by further 
increasing the potential, the low-separation pressures become negative 
again, and finally 4) the isotherm becomes attractive again irrespective of 
the film thickness. The equilibrium contact angle for the low-salinity case 
ranges from 6.40 to 5.33° with subsequent increasing-decreasing trends, 
and for the high-salinity one has an average of 6.02°. 

The results of the sensitivity analysis on the surface charge density are 
also presented in the second section of the supplementary materials file 
for two different charge values. The featured characteristics of this 
analysis are also being discussed. In summary, the equilibrium contact 
angles comprising both surface charge density scenarios range from 4.86 
to 7.83° for the low-salinity case. As can be seen, there is no meaningful 
difference with respect to the contact angle results compared to the 
constant potential case. Thus, from here forward, the constant potential 
solution will be used to present the model results. 


3.3. Impact of temperature, salinity, and interfacial tension 


The equilibrium contact angle is a weak function of the temperature 
and salt concentration (Figs. 4 and 5 in supplementary materials, 
respectively). The equilibrium contact angle demonstrates a decreasing 
trend with increasing temperature from 298 to 380 K. This trend is 
consistent with the experimental observations of Hosseini et al. [22], Ali 
et al. [34], but at odds with the studies of Iglauer et al. [37] (they 
demonstrated an increasing trend), and Hashemi et al. [32] (no trend was 
detected). 

Regarding the impact of salinity, there is a short period of increasing 
trend (concentrations roughly below 0.02 M) followed by a major 
decreasing trend. However, the literature data show opposing salinity 
effects. Similar to the results of this study, the measured CO3 contact 


Table 3 


Upper and lower bounds of the parameters used in the 
sensitivity study. 


Parameter Range 

y, mV —60.00 to 10.00 
o, C/m? —0.006 to 0.006 
n,, (—) 1.40 to 3.00 

P, kPa 20.00 to 130.00 
IFT, mN/m 30.00 to 80.00 
Salinity, M 0.01 to 1.00 


Temperaure, K 298.00 to 380.00 
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Fig. 5. Results of the performed sensitivity to the surface potential at (a) 0.01 M and (b) 1.00 M NaCl concentrations. The varying potential for one of the interfaces 
ranges from —60 to 10 mV and for the other one is equal —15 mV. The other parameters have been reported in the Table 2. 


angles of Ameri et al. [68] demonstrate decreasing salinity effect for 
water-wet surfaces, while the experimental observations of Hosseini et al. 
[22] indicate the direct relationship of contact angle and salinity at 15 
MPa. Hashemi et al. [32] have reported no meaningful trend. It also 
should be mentioned that we have considered only the salinity varia- 
tions, not its influence on other parameters. In real porous media situa- 
tions, changing the salinity affects other parameters, such as surface 
charge density/potential and IFT, and therefore, the trends and the 
impact become different. 

Moreover, this decreasing salinity effect may be inconsistent with the 
LSW observations, where increasing the salinity leads to contact angle 
elevations for water. To further understand the effect of low saline water 
on the contact angle, we have considered the oil-water-rock case study 
presented by Hirasaki [39] and calculated the contact angle for two 
cases: 1) high saline conditions (0.10 M) with high surface potentials 
(—60 mV), and 2) low saline conditions (0.01 M) with lower surface 
potentials (—40 mV). As the salinity decreases, the film will expand and 
from about 0.40 to 6.00 nm, and contact angle decreases from around 15° 
to 7°. It should also be noted that depending on the disjoining pressure 
and capillary pressure curves, there may exist more than one equilibrium 
film thickness as in the case of Hirasaki [39] model, and not all equi- 
librium points reproduce similar salinity effects. 

Returning to the hydrogen problem, the system only has one equi- 
librium point and it corresponds to a decreasing salinity effect as opposed 
to LSW observations. In summary, considering the parameter set and 
assumptions of this study, the effect of salinity on the contact angle can 
also be neglected. Note that this statement does not imply the insignifi- 
cance of salinity in all cases. The salinity effect can be immense when the 
electrostatic forces dominate the vdW ones, for example during LSW. 

The equilibrium contact angle is also an intermediate function of the 
interfacial tension (Fig. 6 in the supplementary materials). The reported 
values of IFT for hydrogen-water in the literature [66,69] range from 42 
to 73 mN/m. Using the Equation (3) and by changing the IFT from 30 to 
80 mN/m, the equilibrium contact angle reduces; in this case, from 8.52 
to 5.21°. The results for the constant charge case are slightly different 
with contact angles falling from 8.39 to 5.13°. 
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Fig. 6. The model response to variations in rock refractive index from 1.40 to 
3.00 at constant potential. The other parameters have been reported in the 
Table 2. (a) represents the vdW disjoining pressure and (b) shows the meniscus 
angle variations. 


3.4. Impact of dielectric constant and refractive index 


The dielectric properties of the system components can also impact 
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the vdW isotherm and the resulting equilibrium and meniscus contact 
angles. The properties of the non-wetting phase (here hydrogen), namely 
relative permittivity and refractive index exhibit minor changes in the 
calculated equilibrium contact angles. The resulting equilibrium contact 
angles are almost equal for different gases, such as CO2, N2, and CH4. On 
the other hand, the low dependency of vdW forces on the gas types can 
also be investigated based on their Hamaker constants. Zeng et al. [70] 
reported the Hamaker values for some gases interacting with pure water 
[70,71]. For comparison, we employed their reported gas Hamaker 
constants and calculated the total one based on the geometric averaging 
rule (Equation (1) in supplementary materials) [39] in a gas-brine-quartz 
system. The results are shown in Table 4. It should be mentioned that the 
rock and water Hamaker values are set to 3.70 x 107”? and 6.50 x 10°” 
J, respectively [39,41]. As it can be seen, the variations are too small to 
have a contributing effect on the vdW force and the resulting contact 
angles. 

This observation may imply that the presence of different gases 
cannot change the nature of vdW interactions. Therefore, as the vdW 
forces dictate the response of the system and because its magnitude is 
almost the same for different gases (Table 4), the final equilibrium con- 
tact angle will not be significantly different for different gases. This 
statement applies as long as the electrostatic interactions remain signif- 
icantly weak compared to the vdW interactions. On this account, one can 
delineate that neutral gases, such as Nz and H3 that can weakly promote 
interactions, would not influence the interaction type and finally, the 
contact angle. However, the presence of polar gaseous components, like 
CO% may promote interactions between the interfaces, and hence, its 
presence may influence the resulting contact angle. 

As mentioned earlier, the second term in the Hamaker constant 
equation (Equation (2) of the supplementary materials) demonstrates the 
London dispersion interactions, which highly dominates the zero- 
frequency term and dictates the interaction type (repulsion/attraction). 
As a result the interactions are mostly influenced by the parameters 
affecting the second term. On this account, the system has very little 
sensitivity towards the rock dielectric constant variations. However, the 
model shows meaningful sensitivity when the refractive index is varied 
(Fig. 6). As can be seen, elevating the refractive index from 1.40 to 3.00 
substantially increases the vdW pressures as well as the meniscus angles. 
Worth noting that the refractive index for quartz and calcite minerals 
usually ranges from 1.40 to 1.90 [41, 72], and as such the considered 
upper range shows an extreme case. The equilibrium contact angle ex- 
hibits 3.00° change, starting from 5.19° and ending at 8.21°. 


4. Effect of roughness on the apparent contact angle 


Since all results presented here are valid for mathematically flat 
surface, they cannot be directly applied to nature systems with rough 
surfaces. To investigate the impact of roughness, Cassie-Baxter theory 
[19] was used. An important point to mention is that the Wenzel equa- 
tion [20] could not be used to study the effect of roughness due to the 
geometry and contact angle limitations (equilibrium contact angles 
driven from the DLVO calculations present low values and hence, their 
incorporation into the Wenzel equation would lead to illogical cosine 


Table 4 

The reported and calculated Hamaker constants for different gases. The Aj; [70] 
represents the Hamaker values of two pure gas layers interacting across a water 
layer, and therefore is always positive (the vdW force between similar media is 
always attractive). The A parameter (Equation (1) in supplementary materials) 
shows a calculated gas-brine-rock Hamaker constant based on the geometric 
averaging rule. 


Gas Type Ay; J [70] A, J ([39]) 

CO 3.15 x 107? — 1.09 x 107” 
N2 3.40 x 107? — 1.08 x 107” 
CH, 3.29 x 107? — 1.09 x 107” 
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values; i.e. > 1). The roughness depth (d) and spacing (w) were taken 
from the study of Sari et al. [31] after post-processing the AFM data 
(Table 5). The roughness ratio (rp) and the wetted fractions (fı and f2) 
were then calculated using the methodology presented in the previous 
section (Fig. 3). 

The influence of the capillary pressure on the contact angle of a 
smooth surface and equilibrium film thickness is shown in Fig. 7 of the 
supplementary materials. The equilibrium film thickness decreases with 
increasing the capillary pressure, and consequently, the equilibrium 
contact angle increases. In Fig. 7, we have assumed a P, — Sy curve from 
the literature for a rock system [35]. Based on the capillary pressure 
value at each saturation, we estimated the radius of curvature and 
apparent contact angle for two sets of protrusion parameters [31]. With 
the rise in the capillary pressure at low saturations, the f parameter for 
the water substantially decreases. This results in considerable increase of 
the corresponding apparent contact angles in water phase. 

For the employed geometry, there is also a strong correlation between 
the apparent contact angle and the roughness parameters, namely the 
spacing and depth. Note that the entry capillary pressure for each pro- 
trusion geometry is different. For instance, increasing d or reducing w 
will increase the entry capillary pressure to the protrusion. It can also be 
seen from Fig. 7b that the geometry with higher w has a smaller entry 
capillary pressure. 

Fig. 8 demonstrates the influence of the groove geometry on the 
apparent contact angles. In this regard, for each d, we have considered a 
w range of around 1-6 pm [31] and calculated the f and rẹ values. It is 
noticeable that rising the spacing at each depth substantially adds to the 
contact angles by above 50°. This is consistent with the observed trend in 
the literature as clearly depicted in the study of Cho and Park [59]. He 
et al. [58] also employed a micro-structure geometry to investigate the 
effect of micro-rod parameters on the wetting state. They concluded that 
depending on the micro-rod height, the wetting can either follow Wenzel, 
Cassie-Baxter, or a transition state. The higher micro-rod heights pre- 
sented the Cassie-Baxter state with contact angles over 90°. For this case, 
increasing the spacing added to the contact angles, consistent with the 
present study. The lower heights, on the contrary, fell into the Wenzel 
state with 0 < 90°. 

Moreover, there exists an inverse correlation between the roughness 
depth and apparent contact angles. The exception would be when the 
roughness depth becomes very large (d = 2 um), where the graph crosses 
the other plots at around w = 3 pm. However, based on the literature 
data, the protrusion depth values mostly peak at below 1 um for rock 
surfaces [22,31] and as such the d = 2 pm plot may consider an extreme 
case. Moreover, it is important to note that the assumed geometry is not 
applicable to all w/d ratios. At some specific w/d ratios, the employed 
geometry cannot accommodate the logical values of fy and therefore, no 
apparent contact angle will be detected by the model. 


5. Validation of the theoretical background 


In this section, a case study will be presented to validate the intro- 
duced methodology for apparent contact angle calculation against the 
available experimental data in the literature. On this account, we make 
use of the data presented by Alnoush et al. [21] for Calcite surface with 
variable surface roughness in presence of oil and water. The RMS 
roughness and apparent contact angles on the rough Calcite surfaces are 
indicated in Table 4 and Fig. 5 of the mentioned study with respect to 


Table 5 
The AFM roughness parameters used in this study. The data are taken from the 
study of [31]. 


RMS Roughness, (nm) Peak-to-Peak Spacing, (nm) Roughness Ratio, (—) 


17.00 435.00 1.00 
336.00 4055.00 1.50 
943.00 5934.00 1.82 
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Fig. 7. The P, — Sy curve [35] and the corresponding apparent contact angles 
calculated in this study. (a) represents a lower roughness geometry with w = 
0.40 um and d = 0.02 um and (b) shows the rougher geometry with w = 4.00 um 
and d = 0.30 pm. The model parameters have also been reported in the Table 2. 
When the P, increases, the contact angle will rise as the liquid will be squeezed 
to the corners. 
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Fig. 8. The apparent contact angle for a hydrogen-brine-rock system as a 
function of protrusion depth and spacing. Increasing the w results in a decrease 
in the total roughness ratio and therefore, the final contact angle will rise. With 
an exception of d = 2 1m, the opposite behavior is seen by changing parameter 
d (refer to Fig. 3). The model parameters have also been reported in the Table 2. 


different diamond grit sizes [21]. The RMS roughness is equal 0.289, 
0.328, and 0.893 um for 1250, 1000, and 600 grit sizes, respectively. The 
naturally cleaved Calcite surface has also a RMS equal 0.030 pm. 

In order to re-generate the apparent contact angles based on the this 
study's methodology, we have utilized the provided RMS roughness and 
adjusted the protrusion spacing at a fixed capillary pressure equal 20 kPa. 
Note that we have used the constant charge solution for electrostatic 
force calculation, and other model parameters include IFT = 25 mN/m, T 
= 293.15 K, qj = —0.004 mN/m, and A = 1 x 107? J [39]. The apparent 
contact angle for the naturally-cleaved calcite surface is equal 70.32° 
with w = 0.17 pm. The results of apparent contact angle calculation for 
the rough Calcite surfaces are also presented in Fig. 9. It can be concluded 
that the presented model can effectively re-produce the experimental 
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Fig. 9. A comparison between apparent contact angles for three Calcite surfaces 
with different roughness (RMS roughness is equal 0.289, 0.328, and 0.893 um 
for 1250, 1000, and 600 grit sizes) calculated in this study and presented by 
Alnoush et al. [21]. The adjusted protrusion spacing for generating the data is 
also presented above each graph. The model parameters are P, = 20 kPa, IFT = 
25 mN/m, T = 293.15 K, q; = —0.004 mN/m, and A = 1 x 107” J [39]. 


observations by adjusting the groove parameters. This is a viable 
approach as different surfaces present unique roughness geometries, and 
even a single substrate may have varying properties at different locations. 


6. Implications for applications 


Contact angle has been vastly used for quantifying wettability of 
underground systems. However, it is been argued that ascribing an 
average contact angle to a whole domain can become troublesome, as 
different pores may present different wetting features. This wetting 
characteristic not only is influenced by the subsurface fluid-rock chem- 
istry but also by the pore geometry and the sub-pore scale surface 
roughness [13,73,74]. Hence, the wetting can be assumed as an 
scale-dependent feature, where its spatial correlation gains importance 
rather than an average value [75-77]. Experimentally-driven contact 
angles on the idealized surfaces cannot also describe the proper wetting 
features at real, small-scale situations. The aim of the current study was 
to understand the physics of the equilibrium contact angle through the 
DLVO theory as well as the apparent contact angle using a 2D 
grooved-shaped roughness. Thus, we employed the Cassie-Baxter theory, 
which could effectively incorporate the required geometries. 

According to the methodology adapted in this paper, increasing the 
roughness in an initially hydrophilic surface tends to decrease the 
apparent contact angle. We investigated the effect of roughness in terms 
of two geometry parameters, namely protrusion spacing (w) and depth 
(d). The overall increase of the roughness was achieved by either 
increasing the d or reducing the w parameters. Particularly, we focused 
on the influence of capillary pressure on the contact angles of both flat 
(Supplementary materials, Fig. 7b) and angular surfaces (Fig. 7). The 
results indicate a positive correlation between these contact angles and 
the capillary pressure. When the capillary increases, the water is 
squeezed to the corners of the groove and as a result, the contact angle 
increases substantially. 

The capillary pressure dependence of contact angle is also been 
investigated in a few studies in literature [13,17,75]. Our observations 
are in line with the numerical study of Armstrong et al. [75], where a 
volume-of-fluid model was utilized for investigating the oil-water-rock 
wettability. Similarly, they observed a noticeable rise and reduction in 
the contact angle and the wetted water fraction. This observation have 
been attributed to the surface aging at higher capillary pressures that 
makes the system less water wet. Riicker et al. [13] also observed the 
influence of global capillary pressure on the surface coverage of the 
non-wetting phase during core-scale experiments. They demonstrated 
that higher capillary pressures lead to higher surface coverage of the oil 
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phase in water-wet rocks. Another important point is that the minimum 
entry capillary pressure is not the same for all geometries. At some 
spacing-to-depth (w/d) ratios, the capillary pressure becomes too high for 
hydrogen to enter the groove. This also highlights that the invasion of a 
pore/roughness does not happen at the same equilibrium contact angle 
(on the flat surface) and it can even present water-repelling features. Sun 
et al. [17] also performed experimental studies and Lattice-Boltzmann 
simulations to investigate air-water wetting at pore scale. They pro- 
vided a link between capillary pressure and contact angle by utilizing the 
P.-S,, data and geometrical models and demonstrated a distribution of 
contact angles at pore scale. 

The results can also imply the impact of phase saturations on the 
observed apparent contact angles [78]. As demonstrated in Fig. 7, in- 
crease of capillary pressures during the drainage process leads to 
reduction of the wetting phase saturation. At high water saturations 
contact angle values do not significant change with the capillary pres- 
sure. As the capillary rises, the apparent contact angles also increase and 
exceed 90°. This can underline the importance of pore geometry and 
roughness near the residual water saturation, where the capillary be- 
comes too high and consequently, the contact angles increase signifi- 
cantly. Accordingly, in the two-phase region (larger than the residual 
values), it may be possible to neglect the influence of the angular ge- 
ometry on the contact angles and on the storage process. 

On the other hand, the fluid-rock interactions is not well-known 
below the residual water saturation. Bear et al. [79] probed into the 
interfacial area in the residual saturation region. They showed that the 
meniscus may represent the shape of a saddle, where the curvature may 
be negative or positive, simultaneously. The change of interface curva- 
ture can also pave the way for the presence of regions with high 
water-repelling apparent contact angles. This observation can be deemed 
consistent with the present Cassie-Baxter model, where the 
water-repelling contact angles represent the transformation of curvature 
from concave to convex. The change in the interface curvature have also 
been seen in the papers of [80-82] using direct numerical simulations, 
although for intermediate-wet regions. 

Finally, according to the Cassie-Baxter theory and our numerical 
procedure, the angular geometry and capillarity can be responsible for 
changing the surface behaviour from hydrophilic (@ < 90°) to hydro- 
phobic (0 > 90°). In our view, the contact angles > 90° may not explicitly 
imply change of the surface wettability towards being gas wet. These may 
present cases where the geometry dictates the distribution of phases 
rather than surface mineralogy/phase types. Consistently, it has been 
seen that some forms of roughness can be used for turning the wetting 
condition of the surface from hydrophilic to water repelling and vice 
versa [24,56—58,83]. 


7. Conclusions 


This study probed the impact of subsurface properties on the wetta- 
bility of the hydrogen-brine-rock system relevant for underground stor- 
age of hydrogen. We utilized the analytical surface force methodology 
and conducted an extensive sensitivity analysis to quantify the impact of 
subsurface parameters on the contact angles. The key contribution of this 
study is the investigation of surface roughness when a water film has 
partially filled the groove. Moreover, we also provided a comparison 
between different gases, like Nz and Hz to examine their differences in 
terms of the storage process. The key conclusions of this study can be 
summarized as: 


@ The constant charge density boundary (for solving the Poisson- 
Boltzmann Equation) resulted in higher equilibrium contact angles 
than the constant potential case. No other difference was identified 
between the results derived from these boundaries. 

@ The total active force and equilibrium contact angles demonstrate 
moderate sensitivity to the interface potentials/charge densities at 


JCIS Open 8 (2022) 100063 


low salinities, and when the salinity rises, the influence of these 
properties on the equilibrium contact angles becomes much weaker. 

@ Considering the sole salinity effect, the equilibrium contact angle 
mainly decreases with increasing salt concentrations [68]. However, 
the literature shows some discrepancies where Hosseini et al. [22] 
reported the direct relation of contact angle and salinity, while 
Hashemi et al. [32] demonstrated no trend. At very high salinity 
conditions, the effect of the electrostatic-driven parameters on the 
wettability can be neglected. 

@ Capillary pressure, as previously reported [13,75], exhibits a positive 
correlation with the contact angles and has an indirect relation with 
the equilibrium film thickness. 

@ Consistent with previous studies [22,32,33], increasing temperature 
and salinity imposed a minuscule decrease to the equilibrium contact 
angles and therefore, their effect can almost be neglected. 

@ The model showed no sensitivity to hydrogen's dielectric properties. 
Rock refractive index imposed large variations to the equilibrium 
contact angle by about 3°. 

@ Contact angle has a strong dependency on surface roughness [58,84, 
85] and capillary pressure [13,75]. The results exhibit a marked in- 
crease in apparent contact angles by rising the protrusion spacing. 
However, the increase in the groove depth results in reduction of 
contact angles, leading the system towards being more water-wet. 

@ During the drainage process, as the capillary increases and interfacial 
area decreases, the water wetted fractions decrease, which result in 
larger apparent contact angles [75]. This indicates the importance of 
roughness and pore geometry only near the region with residual 
water saturation. 

@ The vdW parameters for different gases are similar and hence, contact 
angles cannot be highly affected by these properties. The contact 
angle variations may be achieved only if the water-gas interface 
demonstrates largely-varying properties (for example, in presence of 
polar components, like CO2). 


8. Future work 


Overall, further studies are needed to fully comprehend the influence 
of the surface roughness and pore geometry on the wetting state of the 
subsurface system. For instance, whether severe wetting situations 
observed in this study can also be observed in real subsurface situations. 
The application of the DLVO theory for equilibrium contact angle cal- 
culations can also be modified. It is a question whether hydrogen can 
penetrate the crystal lattice and if its diffusion can impact the electro- 
static properties, which can be captured by other short-range forces. 
Moreover, We have utilized a hypothetical repulsion regarding the non- 
DLVO hydration forces. The presence of different ions with differing 
valencies can also be investigated through a more accurate inclusion of 
hydration forces. 
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Nomenclature 

Letters 

TI Disjoining Pressure 

v Frequency 

€ Dielectric Constant 

p Density 

y Electrical Potential 

K Debye Length 

Aeq Equilibrium Contact Angle 
Oa Apparent Contact Angle 
Om Meniscus Angle 
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Interfacial Tension 
Charge per Electron 
Refractive Index 
Temperature 

Planck's Constant 
Hamaker Constant 
Solution Normality 
Film Thickness 

Ionic Strength 
Charge Density 
Wetted Fraction 
Molar 

Structural Coefficient 
Roughness Depth 
Roughness Frequency 
Boltzmann Constant 
Decay Length 
Capillary Pressure 


rf Roughness Ratio 
Abbreviations 
vdW van der Waals 
HF Hydration Force 
EDL Electrical Double Layer 
PBE Poisson-Boltzmann Eq. 
WH Water-Hydrogen 
RW Rock-Water 
IFT Interfacial Tension 
MW Molecular Weight 
LSW Low Salinity Water 
Subscripts 
r Rock 
w Water 
h Hydrogen 
1 Rock-Water 
2 Water-Hydrogen 
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